library(ggplot2)

rlt_theme <- theme(axis.title.y = element_text(colour="grey20",size=15,face="bold"),
        axis.text.x = element_text(colour="grey20",size=15, face="bold"),
        axis.text.y = element_text(colour="grey20",size=15,face="bold"),  
        axis.title.x = element_text(colour="grey20",size=15,face="bold"))+
        theme_bw()


Beta regression

El método presentado aquí es bastante innovador (2010 en adelante) y no van a encontrar mucha información en la literatura o el web sobre el método. Puede encontrar información suplementaria en el Vignette del paquete betareg.

Referencias:

Cribari-Neto, F., and Zeileis, A. (2010). Beta Regression in R. Journal of Statistical Software, 34(2), 1–24. http://www.jstatsoft.org/v34/i02/.

Grün, B., Kosmidis, I., and Zeileis, A. (2012). Extended Beta Regression in R: Shaken, Stirred, Mixed, and Partitioned. Journal of Statistical Software, 48(11), 1–25. http://www.jstatsoft. org/v48/i11/.

La regresión beta es una acercamiento GLM. La regresión beta es para modelar variables dependientes distribuidas con la distribución beta. Datos que tienen distribución beta incluye proporciones y razones, donde los valores se encuentra entre 0 y 1. Además de producir una regresión que máxima la verosimilitud (tanto para la media como para la precisión de una respuesta distribuida en beta), se proporcionan estimaciones con corrección de sesgo.

Los valores de la variable de respuesta deben estar entre 0 y 1, Si por consecuencia tiene valores que son 0 o 1, es necesario cambiarlo a 1= .999 y 0= 0.001. Los números no pueden ser 0 ni 1, deben ser mayores que 0 y menores que 1.

La paquete “betareg” y sus datos. Tenga en cuenta que este es el enfoque del modelo GLM con la respuesta a través de una función de enlace y un predictor lineal. Igual como GLM normal, existen numerosas funciones de enlace, que pueden ser útil como “logit”, “probit”, “cloglog”, “cauchit”, “log”, “loglog”.

Casi toda la información que se estará presentando aquí proviene de Cribari-Neto y Zeileis (Beta Regression in R).

Consulte el pdf, https://cran.r-project.org/web/packages/betareg/betareg.pdf para obtener una descripción del paquete y más detalles.

Estaré usando datos diferentes a lo que se presenta en el paquete para demostrar el valor de los análisis.

Primer paso, que es una distribución beta?

Lo más importante de la distribución beta es que los valores NUNCA son menores de 0 ni mayores de 1. En adición los intervalos de confianza tampoco no pueden ser ni menor de 0 ni mayor de 1.

library(tidyverse)
library(ggplot2)
library(betareg)

Aquí una serie de distribuciones beta. La distribución beta se calcula con dos parametros, shape 1 o \(\alpha\) y shape 2 o \(\beta\). No vamos a entrar en estos parámetros, aunque pueden ir a la pagina de Wikipedia para tener más información. Lo que deberian fijarse es que aparte de cuando hay un shape que es igual como Shape 1 = 0.5, y shape 2 = 0.5, (tambien shape 2 y 2) todas las otras distribuciones no son simétrica. Siempre hay una cola que se extiende en los valores pequeños o grandes.

library(ggplot2)

ggplot(betadj, aes(b1, cases, colour=beta_shapes))+
  geom_line()+
  ylim(0, 5)+
  scale_colour_discrete(name = "Beta_Shapes", labels = c("Shape 2 y 2", "Shape .5 y .5", "Shape 2 y 5", "Shape 8 y 1.5"))

Propoción de fumadores por diferentes paises. Los datos provienen de World Bank en el siguiente enlace, Smokers. Se encuentra información sobre 187 pais y la proporción de población mayor de 15 años que fuman

library(readr)
Proportion_smokers_world <- read_csv("Data/Proportion_smokers_world.csv")
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   Country_Name = col_character(),
##   Country_Code = col_character(),
##   Indicator_Name = col_character(),
##   Indicator_Code = col_character(),
##   Y2000 = col_double(),
##   Y2005 = col_double(),
##   Y2010 = col_double(),
##   Y2011 = col_double(),
##   Y2012 = col_double(),
##   Y2013 = col_double(),
##   Y2014 = col_double(),
##   Y2015 = col_double(),
##   Y2016 = col_double()
## )
Smokers=Proportion_smokers_world 

Primero vamos a convertir los datos en proporción ya que el programa tiene que utilizar datos mayor de 0 y menor de 1. Seleccionamos el año 2000 y creamos un histograma de la distribución.

Smokers$Y2000P=(Smokers$Y2000)/100 # convertir en proporción

Convertir el promedio y varianza en shape

Convertir el promedio y la varianza de los datos en los valores del shape \(\alpha\) y \(\beta\). Se usa la siguiente para calcular los shapes. El \(\alpha\) referiere al promedio y \(\beta\) a la varianza. Usamos la siguiente función.

estBetaParams <- function(mu, var) {
  alpha <- ((1 - mu) / var - 1 / mu) * mu ^ 2
  beta <- alpha * (1 / mu - 1)
  return(params = list(alpha = alpha, beta = beta))
}
#mean(Smokers$Y2000P)
#var(Smokers$Y2000P)
estBetaParams((mean(Smokers$Y2000P)), (var(Smokers$Y2000P)))
## $alpha
## [1] 3.592488
## 
## $beta
## [1] 9.181559
Smokers$Y2000P=(Smokers$Y2000)/100 # convertir en proporción
x <- seq(0, 1, len = 100)

mean(Smokers$Y2000P)
## [1] 0.2812334
var(Smokers$Y2000P)
## [1] 0.01467551
ggplot(Smokers, aes(Y2000P))+
  geom_histogram(aes(y=..density..), colour="white", fill="grey50")+
  stat_function(aes(x = Smokers$Y2000P, y = ..y..), fun = dbeta, colour="red", n = 100,
                      args = list(shape1 = 3.593, shape2 = 9.185))+
  stat_function(fun = dnorm, 
                args = list(mean = mean(Smokers$Y2000P, na.rm = TRUE), 
                sd = sd(Smokers$Y2000P, na.rm = TRUE)), 
                colour = "green", size = 1)+
  xlab("Proportción de fumadores por Pais")+
  annotate("text", x = .5, y = 4.2, label = "Verde: dist Normal", color="darkgreen")+
  annotate("text", x = .5, y = 3.7, label = "Verde: dist Beta", color="red")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Calculando el intervalo de confianza del promedio de una distribución beta. Se necesita los siguientes paquetes simpleboot , boot.

library(simpleboot)
## Simple Bootstrap Routines (1.1-7)
library(boot) # paquete para calcular el intervalo de confianza de una distribución beta
n=187 # El tamaño de muestra de los datos
alpha = 3.593 # estimado de la función arriba 
beta = 9.185
x = rbeta(n, alpha, beta)


x.boot = one.boot(x, median, R=10^4) # Aquí se usa la mediana, pq el promedio *mean* sera sesgado a la derecha.  
boot.ci(x.boot, type="bca")
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = x.boot, type = "bca")
## 
## Intervals : 
## Level       BCa          
## 95%   ( 0.2385,  0.2893 )  
## Calculations and Intervals on Original Scale

Sobreponemos el intervalo de confianza de la mediana sobre la distribución beta

ggplot(Smokers, aes(Y2000P))+
  geom_histogram(aes(y=..density..), bins=20, colour="white", fill="grey50")+
  stat_function(aes(x = Smokers$Y2000P, y = ..y..), fun = dbeta, colour="red", n = 100,
                      args = list(shape1 = 3.593, shape2 = 9.185))+
  geom_vline(xintercept =0.2320, colour="blue")+ # El intervalo de confianza del promedio
  geom_vline(xintercept =0.2768, colour="blue")+
  rlt_theme+
  xlab("Proporción de fumadores por Pais")

Ahora vamos es hacer el primer análisis de regresión donde nuestra respuesta es una proporción.

library(ggversa) # paquete para los datos
attach(dipodium)
dipodium
Tree NumberTree speciesDBHPlant numberRamet numberDistanceOrientationNumber_of_FlowersHeight_InfloHerbivoryRowPosition_NFNumber_Flowers_positionNumber_of_fruitsPerc_FR_setpardalinum_or_roseumFruit_position_effectFrutos_si_o_noP_or_R_Infl_LenghtNum of fruitsSpecies_NameCardinal orientationPropFR
1E.o75  112.4740  1135n12400     r10r0r10.0001
1E.o76  211.9750  1947n22300     r20r0r20.0001
2E.o76  311.95350  1863n32510.04  r30r1r80.0401
3E.o58  413.24210  2447n42050.25  r40r5r50.25  
4E.o  510.8580  2561n51300     r50r0r20.0001
5E.o59  612.62160  1735n62520.08  p60r2p40.0801
5E.o59  712.82170  1951n72000     p70r0p40.0001
6E.o8  813.12245  1443n82700     r80r0r60.0001
7E.o11.5911.12208  1642n92710.037 r90r1r50.0371
8E.o8.51010.75360  2360n102170.333 p100r1p80.333 
8E.o8.51110.0950  1451n111140.364 p110r1p20.364 
8E.o8.51211.1290  442H122000     r120r0r20.0001
8E.o8.51221.1290  2361n133000     r130r0r20.0001
9E.o12  1311.32135  2150n142850.179 r140r5r30.179 
10E.o73  1415.06210  2153n151320.154 r150r2r50.154 
10E.o73  1425.06210  1635n162400     r160r0r50.0001
10E.o73  1435.06210  H172450.208 r170r5r50.208 
11E.o73  1514.16322  1344n182810.0357p180r1p80.0358
12E.o16  1611.32125  3379n192640.154 r190r4r30.154 
12E.o16  1611.32125  1560n20900     r200r0r30.0001
13E.o41  1710.52120  1743n211860.333 r210r6r30.333 
14E.o14  1812.5630  1443n221920.105 p220r2r10.105 
14E.o14  1911.05330  2352n231920.105 r230r2p80.105 
15E.o11  2014.9 165  3069n243030.1   r240r3r40.1   
15E.o11  2014.22170  2862n252920.069 p10r2p40.0691
16E.o84  2111.66220  1845n261900     r20r0r50.0001
17E.o79  2216.5 130  1750n272200     r30r0r30.0001
17E.o79  2313   115  3563n282050.25  r40r5r30.25  
17E.o79  2412.8 115  1334n293230.0938r50r3r30.0939
17E.o79  2511.3 215  3073n301320.154 r60r2r50.154 
18E.o79.52610.5 250  3267n312730.111 r70r3r60.111 
19E.o84  2712.2220  2165n321130.273 r80r3r10.273 
20E.o84  2812.3880  2872n333200     r90r0r20.0001
20E.o84  2913.0380  1446n342330.13  r100r3r20.131 
21E.radiata7  3011.3 254  1451n351330.231 r110r3r60.231 
21E.radiata7  3113.130  1641n361620.125 r120r2r10.125 
22E.o40  3216.15340  2448n372210.0455r130r1r80.0456
22E.o40  3310.34140  1545n382850.179 r140r5r40.179 
22E.o40  3414.6 110  2649n392010.05  r150r1r30.0501
23E.o80  3515.2150  2650n401500     r160r0p20.0001
23E.o80  3612.45310  2047n411300     r170r0r70.0001
23E.o80  3712.48290  1446n421500     p180r0r70.0001
24E.o68  3811.78140  2348n431910.0526r190r1r40.0527
25E.o69  392.34320  2457n442850.179 r200r5r80.179 
25E.o69  3922.34320  2156n451850.278 r210r5r80.278 
25E.o69  4012.08350  1941n462150.238 r220r5r80.238 
25E.o69  4111.980  2652n472830.107 r230r3r10.107 
25E.o69  4212.645  2243n482780.296 r11r8r10.296 
25E.o69  4312.8115  1935n491120.182 r20r2r10.182 
26E.o22  443.5160  337H501200     r30r0r20.0001
26E.o22  4423.5160  3672n511300     r40r0r20.0001
27E.o21.54513.7 21.53250n523530.0857r50r3r10.0858
28E.o39  462.9250  2152n532500     r60r0r20.0001
28E.o39  4622.9250  2351n542660.231 r70r6r20.231 
28E.o39  4712.46352  2352n551710.0588r80r1r80.0589
28E.o39  4812.46240  1639n562620.0769r90r2r60.077 
29E.o14  4910.6810  2263n572210.0455r100r1r10.0456
29E.o14  5011.21300  1342n582100     r110r0r70.0001
29E.o14  5111.11245  2448n592700     r120r0r60.0001
29E.o14  5211.25155  2865n601600     r130r0r40.0001
30E.o15  5311.6 290  2871n611740.235 r140r0r70.235 
31E.o83  5413   300  2454n623220.0625r150r1r70.0626
31E.o83  5512.71270  531H631600     r160r2r60.0001
32E.o44  5611.9 10  930n     170r0r1     
32E.o44  5711.7910  3777n     180r1     
32E.o44  5811.525  2460n     190r1     
32E.o44  590.95310  1859n     200r7     
32E.o44  5920.95310  2761n     210r7     
32E.o44  6011   310  2971n     220r7     
33E.del57  6113.14290  2455n     230r7     
34E.o79  6213.67190  2439n     240r5     
35E.o10  6313.7230  2962n     11r1     
35E.o10  6414.88125  1648n     20r3     
36E.o79  6512.8 295  2452n     31r7     
37E.o8.56610.69142  18n     41r4     
37E.o8.56710.44132  5n     51r3     
38E.o84  6811.31140  2062n     60r4     
38E.o84  6921.31140  1230n     71r4     
39E.o41  7011.86240  1943n     80r6     
41E.o3.57111.36120  2555n     90r3     
40E.o13  7211.2 240  n     100r6     
40E.o13  7221.2 240  n     110r6     
40E.o13  7231.2 240  n     120r6     
42E.del7.57311.25120  3280n     130r3     
42E.del7.57321.25120  2157n     140r3     
42E.del7.57331.25120  2656n     150r3     
42E.del7.57341.25120  2964n     160r3     
42E.del7.57351.25120  2262n     170r3     
42E.del7.57361.25120  1444n     180r3     
42E.del7.57411.53235  2253n     190r6     
42E.del7.57511.66240  2150n     200r6     
42E.del7.57611.61245  2242n     250r6     
42E.del7.57621.61245  2152n     10r6     
42E.del7.57631.61245  841n     20r6     
42E.del7.57711.3 210  1133n     30r5     
42E.del7.57721.3 210  1746n     40r5     
42E.del7.57731.3 210  2468n     50r5     
42E.del7.57811.64210  1448n     60r5     
42E.del7.57821.64210  2361n     70r5     
42E.del7.57831.64210  2350n     80r5     
42E.del7.57911.72210  1747n     90r5     
42E.del7.57921.72210  1039n     100r5     
42E.del7.58011.78180  3377n     110r4     
42E.del7.58112.19180  1440n     120r4     
42E.del7.58212.32180  1235n     130r4     
42E.del7.58222.32180  2538n     10p4     
42E.del7.58312.42180  2871n     20p4     
42E.del7.58322.42180  2766n     30p4     
42E.del7.58412.55190  2347n     40p5     
42E.del7.58512.61174  1229n     50p4     
42E.del7.58612.87174  1762n     60p4     
42E.del7.58712.63170  1851n     70p4     
42E.del7.58812.39170  2345n     80p4     
42E.del7.58911.54158  2251n     90p4     
42E.del7.59012.14152  2558n     100p4     
42E.del7.59111.95150  937n     110p4     
42E.del7.59211.43136  1853n     120p4     
42E.del7.59221.43136  H     131p4     
42E.del7.59231.43136  H     140p4     
42E.del7.59241.43136  H     150p4     
42E.del7.59313.12145  2561n     160p4     
42E.del7.59323.12145  1252n     171p4     
43E.del7.59412.47228  1317n     180p6     
43E.del7.59422.47228  H     190p6     
43E.del7.59432.47228  H     200p6     
43E.del7.59512.76228  H     210p6     
43E.del7.59613.72216  2771n     220p5     
43E.del7.59623.72216  H     230p5     
44E.o64  9713.48340  1130n     240p8     
44E.o64  9812.49350  1746n     250p8     
44E.o64  9822.49350  H     10p8     
45E.o20.59911.36280  3278n     20p7     
45E.o20.510011.12250  1536n     30p6     
45E.o20.510021.12250  2560n     40p6     
45E.o20.510031.12250  2761n     50p6     
45E.o20.510041.12250  2360n     60p6     
45E.o20.510111.1470  2956n     70p2     
46E.o19.510214.18185  2452n     80p5     
46E.o19.510224.18185  H     90p5     
46E.o19.510313.33180  4372n     100p4     
46E.o19.510413.23180  2744n     110p4     
46E.o19.510512.62178  3360n     120p4     
46E.o19.510613.06174  2756n     130p4     
46E.o19.510713.3 167  2564n     140p4     
46E.o19.510813.34164  1640n     150p4     
46E.o19.510913.22160  2968n     160p4     
46E.o19.510923.22160  H     170p4     
46E.o19.511013.53160  1738n     180p4     
46E.o19.511023.53160  3161n     190p4     
46E.o19.511113.48158  1237n     200p4     
46E.o19.511123.48158  726n     10r4     
46E.o19.511211.24160  1748n     20r4     
46E.o19.511221.24160  3578n     30r4     
46E.o19.511311   140  1946n     40r4     
46E.o19.511410.96108  3263n     50r3     
46E.o19.511511.83107  1841n     60r3     
46E.o19.511521.83107  2251n     70r3     
46E.o19.511611.18104  1240n     80r3     
46E.o19.511710.66330  1019n     90r8     
46E.o19.511720.66330  1531n     100r8     
46E.o19.511730.66330  1636n     110r8     
46E.o19.511740.66330  H     120r8     
47E.o62  11813.915  2151n     130r1     
47E.o62  11913.858  1445n     140r1     
47E.o62  12023.858  1644n     150r1     
47E.o62  12033.858  1228n     160r1     
47E.o62  12014.1720  3477n     170r1     
47E.o62  12113.2 350  1434n     180r8     
47E.o62  12213.350  1740n     190r1     
47E.o62  12313.1 324  2865n     200r8     
47E.o62  12413.08324  2032n     210r8     
47E.o62  12513.4 324  2359n     220r8     
47E.o62  12523.4 324  1348n     230r8     
47E.o62  12533.4 324  1646n     240r8     
47E.o62  12613   332  2733n     250r8     
47E.o62  12711.9 305  1133n     260r7     
47E.o62  12811.65290  1537n     270r7     
47E.o62  12911.88290  2867n     10r7     
47E.o62  13012.45284  17n     20r7     
47E.o62  13112.46290  1134n     30r7     
47E.o62  13213.53284  635n     40r7     
47E.o62  13223.53284  1442n     50r7     
47E.o62  13233.53284  2866n     60r7     
47E.o62  13243.53284  2663n     70r7     
47E.o62  13313.53284  3062n     80r7     
47E.o62  13411.11285  1852n     90r7     
47E.o62  13511.57250  2668n     100r6     
47E.o62  13611.55220  3371n     110r5     
47E.o62  13712.62220  2050n     121r5     
47E.o62  13722.62220  2057n     130r5     
47E.o62  13732.62220  2358n     140r5     
47E.o62  13742.62220  621n     150r5     
47E.o62  13752.62220  3065n     160r5     
47E.o62  13812.4 170  2053n     170r4     
47E.o62  13912.65166  2150n     180r4     
47E.o62  14011.75146  1950n     190r4     
47E.o62  14113   143  1947n     200r4     
47E.o62  14213.4 143  3170n     210r4     
            220r     
            230r     
            240r     
            250r     
            260r     
            270r     
            10p     
            20p     
            30p     
            40p     
            50p     
            60p     
            71p     
            80p     
            90p     
            100p     
            110p     
            120p     
            130p     
            140p     
            150p     
            160p     
            170p     
            180p     
            190p     
            200p     
            210p     
            10p     
            20p     
            30p     
            41p     
            50p     
            60p     
            70p     
            80p     
            90p     
            100p     
            110p     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            170r     
            180r     
            190r     
            200r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            170r     
            180r     
            190r     
            200r     
            210r     
            220r     
            230r     
            240r     
            250r     
            260r     
            270r     
            280r     
            290r     
            300r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            121r     
            130r     
            140r     
            150r     
            160r     
            170r     
            181r     
            190r     
            200r     
            210r     
            221r     
            230r     
            241r     
            251r     
            260r     
            270r     
            280r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            101r     
            111r     
            120r     
            130r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            170r     
            180r     
            190r     
            200r     
            210r     
            220r     
            230r     
            240r     
            10r     
            21r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            91r     
            101r     
            111r     
            120r     
            130r     
            140r     
            150r     
            160r     
            171r     
            180r     
            190r     
            200r     
            210r     
            220r     
            230r     
            240r     
            10p     
            20p     
            30p     
            40p     
            50p     
            60p     
            70p     
            80p     
            90p     
            100p     
            110p     
            120p     
            130p     
            140p     
            150p     
            160p     
            170p     
            180p     
            190p     
            200p     
            210p     
            220p     
            231p     
            240p     
            250p     
            260p     
            270p     
            280p     
            11r     
            21r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            170r     
            180r     
            190r     
            200r     
            210r     
            221r     
            230r     
            241r     
            250r     
            260r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            11r     
            21r     
            31r     
            41r     
            51r     
            61r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            170r     
            180r     
            90r     
            10r     
            20r     
            31r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            171r     
            180r     
            190r     
            10p     
            20p     
            30p     
            40p     
            50p     
            60p     
            70p     
            80p     
            91p     
            101p     
            110p     
            120p     
            130p     
            140p     
            150p     
            160p     
            170p     
            180p     
            190p     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            170r     
            180r     
            191r     
            200r     
            210r     
            220r     
            230r     
            240r     
            251r     
            261r     
            270r     
            280r     
            290r     
            300r     
            10p     
            20p     
            30p     
            40p     
            50p     
            60p     
            71p     
            80p     
            90p     
            100p     
            110p     
            120p     
            130p     
            141p     
            150p     
            160p     
            170p     
            180p     
            190p     
            200p     
            210p     
            220p     
            230p     
            240p     
            250p     
            260p     
            270p     
            280p     
            290p     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            170r     
            180r     
            190r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            170r     
            180r     
            190r     
            200r     
            210r     
            220r     
            11r     
            20r     
            31r     
            40r     
            50r     
            61r     
            71r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            171r     
            180r     
            190r     
            200r     
            11r     
            20r     
            31r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            170r     
            180r     
            190r     
            200r     
            210r     
            220r     
            231r     
            240r     
            250r     
            260r     
            270r     
            280r     
            290r     
            300r     
            310r     
            320r     
            11r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            81r     
            90r     
            100r     
            110r     
            120r     
            130r     
            10r     
            21r     
            30r     
            41r     
            51r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            170r     
            180r     
            190r     
            200r     
            210r     
            220r     
            230r     
            240r     
            250r     
            260r     
            270r     
            11r     
            20r     
            31r     
            41r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            170r     
            180r     
            190r     
            200r     
            210r     
            220r     
            230r     
            240r     
            250r     
            260r     
            270r     
            280r     
            290r     
            300r     
            310r     
            320r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            131r     
            140r     
            150r     
            160r     
            170r     
            180r     
            190r     
            200r     
            210r     
            221r     
            231r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            71r     
            80r     
            90r     
            100r     
            111r     
            121r     
            130r     
            10r     
            21r     
            30r     
            40r     
            51r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            81r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            170r     
            180r     
            190r     
            200r     
            210r     
            220r     
            10r     
            20r     
            30r     
            40r     
            50r     
            61r     
            70r     
            80r     
            90r     
            100r     
            111r     
            120r     
            131r     
            140r     
            151r     
            160r     
            170r     
            180r     
            190r     
            200r     
            210r     
            220r     
            230r     
            240r     
            251r     
            260r     
            270r     
            280r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            81r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            170r     
            180r     
            190r     
            200r     
            10p     
            20p     
            30p     
            40p     
            50p     
            60p     
            70p     
            80p     
            90p     
            100p     
            110p     
            120p     
            130p     
            140p     
            150p     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            10r     
            20r     
            30r     
            40r     
            51r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            170r     
            180r     
            190r     
            10r     
            21r     
            30r     
            40r     
            50r     
            60r     
            71r     
            80r     
            91r     
            100r     
            110r     
            121r     
            131r     
            140r     
            150r     
            160r     
            170r     
            180r     
            190r     
            200r     
            210r     
            220r     
            230r     
            240r     
            250r     
            260r     
            270r     
            280r     
            11r     
            21r     
            30r     
            41r     
            50r     
            60r     
            70r     
            81r     
            91r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            170r     
            180r     
            10r     
            21r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            151r     
            160r     
            170r     
            181r     
            191r     
            201r     
            210r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            81r     
            90r     
            100r     
            111r     
            120r     
            130r     
            141r     
            150r     
            160r     
            170r     
            180r     
            190r     
            200r     
            210r     
            220r     
            230r     
            240r     
            250r     
            260r     
            270r     
            280r     
            10r     
            20r     
            30r     
            41r     
            51r     
            60r     
            71r     
            81r     
            90r     
            100r     
            110r     
            121r     
            131r     
            141r     
            151r     
            160r     
            170r     
            180r     
            190r     
            200r     
            210r     
            220r     
            230r     
            240r     
            250r     
            260r     
            270r     
            10r     
            21r     
            30r     
            41r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            161r     
            170r     
            180r     
            190r     
            201r     
            210r     
            220r     
            231r     
            240r     
            250r     
            260r     
            270r     
            280r     
            290r     
            300r     
            310r     
            320r     
            330r     
            340r     
            350r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            170r     
            180r     
            190r     
            200r     
            210r     
            220r     
            230r     
            240r     
            250r     
            10r     
            21r     
            30r     
            41r     
            51r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            121r     
            130r     
            141r     
            150r     
            160r     
            170r     
            180r     
            191r     
            200r     
            210r     
            220r     
            230r     
            240r     
            250r     
            260r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            111r     
            120r     
            130r     
            140r     
            150r     
            160r     
            170r     
            10r     
            20r     
            30r     
            40r     
            51r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            131r     
            140r     
            150r     
            160r     
            170r     
            180r     
            190r     
            200r     
            210r     
            220r     
            230r     
            240r     
            250r     
            260r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            91r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            170r     
            180r     
            190r     
            200r     
            210r     
            220r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            170r     
            180r     
            190r     
            200r     
            210r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            170r     
            180r     
            190r     
            200r     
            210r     
            220r     
            230r     
            240r     
            250r     
            260r     
            270r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            10r     
            20r     
            30r     
            41r     
            50r     
            60r     
            70r     
            80r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            101r     
            110r     
            120r     
            130r     
            140r     
            151r     
            160r     
            170r     
            180r     
            190r     
            200r     
            210r     
            220r     
            230r     
            240r     
            250r     
            260r     
            270r     
            280r     
            290r     
            300r     
            310r     
            320r     
            10r     
            20r     
            30r     
            40r     
            50r     
            60r     
            70r     
            80r     
            90r     
            100r     
            110r     
            120r     
            130r     
            140r     
            150r     
            160r     
library(betareg) # El paquete para hacer regresión beta

Regresión beta

Los datos provienen de una especie de orquidea de Autralia, Dipodium roseum, recolectado por RLT en 2004-2005. Vamos a evaluar la relación entre el número de flores y la proporción de frutos por planta. El primer paso es asegarar que no haya valores de 0 y 1. En este caso no hay ni una planta que tiene 100% de frutos, pero si hay individuos que tienen cero frutos.

  • RECUERDA x >0 y <1. NO se acepta 0 o 1. Entonces a los valores de 0 se le puede añadir un valor mínimo como 0.001 y a los 1 restar 0.001. En realidad esta modificación no impacta la interpretación de los resultados.

  • se remueva también del archivo los NA.

  • Nota que el modelo se construye como un modelo lineal betareg(y~x, data =na.omit(df)). Las variables del archivo son PropFR, el proporción de frutos (número de frutos/números de flores) por cada individuo y el número de flores, Number_of_Flowers.

head(dipodium)
Tree NumberTree speciesDBHPlant numberRamet numberDistanceOrientationNumber_of_FlowersHeight_InfloHerbivoryRowPosition_NFNumber_Flowers_positionNumber_of_fruitsPerc_FR_setpardalinum_or_roseumFruit_position_effectFrutos_si_o_noP_or_R_Infl_LenghtNum of fruitsSpecies_NameCardinal orientationPropFR
1E.o75112.47401135n12400   r10r0r10.0001
1E.o76211.97501947n22300   r20r0r20.0001
2E.o76311.953501863n32510.04r30r1r80.0401
3E.o58413.242102447n42050.25r40r5r50.25  
4E.o510.85802561n51300   r50r0r20.0001
5E.o59612.621601735n62520.08p60r2p40.0801
library(tidyverse)
dipodium$PropFR=dipodium$Perc_FR_set+0.0001 # solucionar para remover los cero
#dipodium$PropFR
dipodium2=dipodium %>% 
  select(PropFR, Number_of_Flowers,Height_Inflo, Distance) %>% 
  filter(complete.cases(PropFR, Number_of_Flowers,Height_Inflo, Distance))
#dipodium2
write.csv(x=dipodium2, file="dipodium2.csv")


library(readr)
dipodium2 <- read_csv("Data/dipodium2.csv")
## Warning: Missing column names filled in: 'X1' [1]
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   X1 = col_double(),
##   PropFR = col_double(),
##   Number_of_Flowers = col_double(),
##   Height_Inflo = col_double(),
##   Distance = col_double()
## )
modelpropFr=betareg(PropFR~Number_of_Flowers+Height_Inflo+Distance, data =na.omit(dipodium2))
summary(modelpropFr)
## 
## Call:
## betareg(formula = PropFR ~ Number_of_Flowers + Height_Inflo + Distance, 
##     data = na.omit(dipodium2))
## 
## Standardized weighted residuals 2:
##     Min      1Q  Median      3Q     Max 
## -3.8613 -0.4046  0.2672  0.8535  1.2798 
## 
## Coefficients (mean model with logit link):
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -3.796016   0.727139  -5.220 1.78e-07 ***
## Number_of_Flowers  0.076200   0.028247   2.698  0.00698 ** 
## Height_Inflo       0.006054   0.017847   0.339  0.73444    
## Distance          -0.111010   0.093825  -1.183  0.23674    
## 
## Phi coefficients (precision model with identity link):
##       Estimate Std. Error z value Pr(>|z|)    
## (phi)   4.8382     0.9982   4.847 1.25e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Type of estimator: ML (maximum likelihood)
## Log-likelihood: 113.9 on 5 Df
## Pseudo R-squared: 0.2757
## Number of iterations: 19 (BFGS) + 2 (Fisher scoring)

Visualizar un grafico de regresión beta

dipodiumbeta=dipodium2[,c("Number_of_Flowers","PropFR")] # crear un df con solamente las columnas de interes 

 dp2=dipodiumbeta[complete.cases(dipodiumbeta),] # remover los "NA" 
 
modelpropFr=betareg(PropFR~Number_of_Flowers, data=dp2) # El modelo con solamente la explicativa
summary(modelpropFr)
## 
## Call:
## betareg(formula = PropFR ~ Number_of_Flowers, data = dp2)
## 
## Standardized weighted residuals 2:
##     Min      1Q  Median      3Q     Max 
## -3.3221 -0.3181  0.2414  0.8225  1.2346 
## 
## Coefficients (mean model with logit link):
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -3.88360    0.44251  -8.776  < 2e-16 ***
## Number_of_Flowers  0.08259    0.01803   4.581 4.62e-06 ***
## 
## Phi coefficients (precision model with identity link):
##       Estimate Std. Error z value Pr(>|z|)    
## (phi)   4.6947     0.9699    4.84  1.3e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Type of estimator: ML (maximum likelihood)
## Log-likelihood: 112.9 on 3 Df
## Pseudo R-squared: 0.2538
## Number of iterations: 18 (BFGS) + 1 (Fisher scoring)
 #head(dp2)
#dp2$PropFR
predict(modelpropFr, type = "response") # calcular los valores estimados (predichos)
##          1          2          3          4          5          6          7 
## 0.02568523 0.02783519 0.03015954 0.04856508 0.05679268 0.05679268 0.05679268 
##          8          9         10         11         12         13         14 
## 0.06138242 0.06138242 0.06138242 0.06138242 0.06138242 0.06138242 0.06631700 
##         15         16         17         18         19         20         21 
## 0.06631700 0.07161801 0.07161801 0.07161801 0.07161801 0.07730766 0.07730766 
##         22         23         24         25         26         27         28 
## 0.07730766 0.08340872 0.08340872 0.08994433 0.08994433 0.08994433 0.08994433 
##         29         30         31         32         33         34         35 
## 0.09693789 0.10441284 0.10441284 0.10441284 0.10441284 0.10441284 0.11239245 
##         36         37         38         39         40         41         42 
## 0.11239245 0.12089957 0.12089957 0.12089957 0.12089957 0.12089957 0.12089957 
##         43         44         45         46         47         48         49 
## 0.12995632 0.12995632 0.12995632 0.12995632 0.12995632 0.13958381 0.14980172 
##         50         51         52         53         54         55         56 
## 0.14980172 0.14980172 0.17207828 0.17207828 0.17207828 0.17207828 0.19690033 
##         57         58         59         60         61         62 
## 0.19690033 0.22433274 0.22433274 0.23903102 0.27035714 0.28695536
dp2$response=predict(modelpropFr, type = "response")
#dp2$link=predict(modelpropFr, type = "link")
dp2$precision=predict(modelpropFr, type = "precision")
dp2$variance=predict(modelpropFr, type = "variance")
dp2$quantile_.01=predict(modelpropFr, type = "quantile", at = c(0.01))
dp2$quantile_.05=predict(modelpropFr, type = "quantile", at = c(0.05))
dp2$quantile_.10=predict(modelpropFr, type = "quantile", at = c(0.10))
dp2$quantile_.15=predict(modelpropFr, type = "quantile", at = c(0.15))
dp2$quantile_.20=predict(modelpropFr, type = "quantile", at = c(0.20))
dp2$quantile_.25=predict(modelpropFr, type = "quantile", at = c(0.25))
dp2$quantile_.30=predict(modelpropFr, type = "quantile", at = c(0.30))
dp2$quantile_.35=predict(modelpropFr, type = "quantile", at = c(0.35))
dp2$quantile_.40=predict(modelpropFr, type = "quantile", at = c(0.40))
dp2$quantile_.45=predict(modelpropFr, type = "quantile", at = c(0.45))
dp2$quantile_.50=predict(modelpropFr, type = "quantile", at = c(0.50))
dp2$quantile_.55=predict(modelpropFr, type = "quantile", at = c(0.55))
dp2$quantile_.60=predict(modelpropFr, type = "quantile", at = c(0.60))
dp2$quantile_.65=predict(modelpropFr, type = "quantile", at = c(0.65))
dp2$quantile_.70=predict(modelpropFr, type = "quantile", at = c(0.70))
dp2$quantile_.75=predict(modelpropFr, type = "quantile", at = c(0.75))
dp2$quantile_.80=predict(modelpropFr, type = "quantile", at = c(0.80))
dp2$quantile_.85=predict(modelpropFr, type = "quantile", at = c(0.85))
dp2$quantile_.90=predict(modelpropFr, type = "quantile", at = c(0.90))
dp2$quantile_.95=predict(modelpropFr, type = "quantile", at = c(0.95))
dp2$quantile_.99=predict(modelpropFr, type = "quantile", at = c(0.99))
dp2
Number_of_FlowersPropFRresponseprecisionvariancequantile_.01quantile_.05quantile_.10quantile_.15quantile_.20quantile_.25quantile_.30quantile_.35quantile_.40quantile_.45quantile_.50quantile_.55quantile_.60quantile_.65quantile_.70quantile_.75quantile_.80quantile_.85quantile_.90quantile_.95quantile_.99
30.00010.02574.690.004393.86e-182.42e-127.58e-102.19e-082.38e-071.51e-066.86e-062.46e-057.46e-050.0001980.0004750.001050.002170.004240.007920.01430.02530.04440.07920.1530.342
40.00010.02784.690.004757.44e-171.66e-113.34e-097.44e-086.72e-073.71e-061.5e-05 4.87e-050.0001350.0003330.0007480.001550.003040.005650.0101 0.01750.02990.05060.08720.1630.353
50.00010.03024.690.005141.14e-159.81e-111.31e-082.3e-07 1.75e-068.48e-063.07e-059.13e-050.0002350.0005390.00114 0.002240.004160.0074 0.0127 0.02120.03490.05720.09570.1730.364
110.00010.04864.690.008112.75e-103.2e-07 6.69e-063.96e-050.00014 0.0003720.0008290.00163 0.00295 0.00497 0.00795 0.0122 0.0182 0.0265 0.0378 0.05340.07510.106 0.155 0.2420.434
130.03580.05684.690.009415.29e-092.21e-062.98e-050.0001360.0004010.0009280.00184 0.0033  0.00547 0.00859 0.0129  0.0187 0.0265 0.0368 0.0503 0.06830.09260.127 0.177 0.2670.458
130.09390.05684.690.009415.29e-092.21e-062.98e-050.0001360.0004010.0009280.00184 0.0033  0.00547 0.00859 0.0129  0.0187 0.0265 0.0368 0.0503 0.06830.09260.127 0.177 0.2670.458
130.00010.05684.690.009415.29e-092.21e-062.98e-050.0001360.0004010.0009280.00184 0.0033  0.00547 0.00859 0.0129  0.0187 0.0265 0.0368 0.0503 0.06830.09260.127 0.177 0.2670.458
140.00010.06144.690.0101 1.96e-085.21e-065.77e-050.0002360.0006410.00139 0.00263 0.00451 0.00723 0.011   0.016   0.0227 0.0314 0.0427 0.0574 0.07650.102 0.137 0.189 0.28 0.47 
140.364 0.06144.690.0101 1.96e-085.21e-065.77e-050.0002360.0006410.00139 0.00263 0.00451 0.00723 0.011   0.016   0.0227 0.0314 0.0427 0.0574 0.07650.102 0.137 0.189 0.28 0.47 
140.105 0.06144.690.0101 1.96e-085.21e-065.77e-050.0002360.0006410.00139 0.00263 0.00451 0.00723 0.011   0.016   0.0227 0.0314 0.0427 0.0574 0.07650.102 0.137 0.189 0.28 0.47 
140.131 0.06144.690.0101 1.96e-085.21e-065.77e-050.0002360.0006410.00139 0.00263 0.00451 0.00723 0.011   0.016   0.0227 0.0314 0.0427 0.0574 0.07650.102 0.137 0.189 0.28 0.47 
140.231 0.06144.690.0101 1.96e-085.21e-065.77e-050.0002360.0006410.00139 0.00263 0.00451 0.00723 0.011   0.016   0.0227 0.0314 0.0427 0.0574 0.07650.102 0.137 0.189 0.28 0.47 
140.00010.06144.690.0101 1.96e-085.21e-065.77e-050.0002360.0006410.00139 0.00263 0.00451 0.00723 0.011   0.016   0.0227 0.0314 0.0427 0.0574 0.07650.102 0.137 0.189 0.28 0.47 
150.00010.06634.690.0109 6.54e-081.15e-050.0001070.0003920.0009890.00203 0.00366 0.00605 0.00937 0.0138  0.0197  0.0273 0.037  0.0493 0.0651 0.08540.112 0.149 0.202 0.2930.483
150.179 0.06634.690.0109 6.54e-081.15e-050.0001070.0003920.0009890.00203 0.00366 0.00605 0.00937 0.0138  0.0197  0.0273 0.037  0.0493 0.0651 0.08540.112 0.149 0.202 0.2930.483
160.03710.07164.690.0117 1.99e-072.39e-050.0001880.0006270.00148 0.00288 0.00499 0.00794 0.0119  0.0172  0.0239  0.0324 0.0431 0.0565 0.0733 0.09480.123 0.16  0.215 0.3070.496
160.00010.07164.690.0117 1.99e-072.39e-050.0001880.0006270.00148 0.00288 0.00499 0.00794 0.0119  0.0172  0.0239  0.0324 0.0431 0.0565 0.0733 0.09480.123 0.16  0.215 0.3070.496
160.125 0.07164.690.0117 1.99e-072.39e-050.0001880.0006270.00148 0.00288 0.00499 0.00794 0.0119  0.0172  0.0239  0.0324 0.0431 0.0565 0.0733 0.09480.123 0.16  0.215 0.3070.496
160.077 0.07164.690.0117 1.99e-072.39e-050.0001880.0006270.00148 0.00288 0.00499 0.00794 0.0119  0.0172  0.0239  0.0324 0.0431 0.0565 0.0733 0.09480.123 0.16  0.215 0.3070.496
170.08010.07734.690.0125 5.56e-074.69e-050.0003170.00097 0.00215 0.00399 0.00664 0.0102  0.015   0.021   0.0286  0.0381 0.0498 0.0642 0.0822 0.105 0.134 0.172 0.228 0.3210.509
170.333 0.07734.690.0125 5.56e-074.69e-050.0003170.00097 0.00215 0.00399 0.00664 0.0102  0.015   0.021   0.0286  0.0381 0.0498 0.0642 0.0822 0.105 0.134 0.172 0.228 0.3210.509
170.00010.07734.690.0125 5.56e-074.69e-050.0003170.00097 0.00215 0.00399 0.00664 0.0102  0.015   0.021   0.0286  0.0381 0.0498 0.0642 0.0822 0.105 0.134 0.172 0.228 0.3210.509
180.04010.08344.690.0134 1.44e-068.76e-050.0005150.00145 0.00304 0.00541 0.00868 0.013   0.0185  0.0254  0.034   0.0444 0.0571 0.0726 0.0916 0.115 0.145 0.185 0.242 0.3350.522
180.00010.08344.690.0134 1.44e-068.76e-050.0005150.00145 0.00304 0.00541 0.00868 0.013   0.0185  0.0254  0.034   0.0444 0.0571 0.0726 0.0916 0.115 0.145 0.185 0.242 0.3350.522
190.00010.08994.690.0144 3.45e-060.0001560.0008080.00212 0.0042  0.00718 0.0112  0.0163  0.0226  0.0305  0.0399  0.0514 0.0651 0.0816 0.102  0.126 0.158 0.198 0.256 0.35 0.535
190.00010.08994.690.0144 3.45e-060.0001560.0008080.00212 0.0042  0.00718 0.0112  0.0163  0.0226  0.0305  0.0399  0.0514 0.0651 0.0816 0.102  0.126 0.158 0.198 0.256 0.35 0.535
190.238 0.08994.690.0144 3.45e-060.0001560.0008080.00212 0.0042  0.00718 0.0112  0.0163  0.0226  0.0305  0.0399  0.0514 0.0651 0.0816 0.102  0.126 0.158 0.198 0.256 0.35 0.535
190.182 0.08994.690.0144 3.45e-060.0001560.0008080.00212 0.0042  0.00718 0.0112  0.0163  0.0226  0.0305  0.0399  0.0514 0.0651 0.0816 0.102  0.126 0.158 0.198 0.256 0.35 0.535
200.00010.09694.690.0154 7.76e-060.0002670.00123 0.003   0.00568 0.00935 0.0141  0.0201  0.0273  0.0361  0.0465  0.059  0.0737 0.0913 0.112  0.138 0.17  0.212 0.27  0.3650.549
210.179 0.104 4.690.0164 1.64e-050.0004380.00181 0.00415 0.00752 0.012   0.0176  0.0244  0.0326  0.0424  0.0538  0.0673 0.0831 0.102  0.124  0.151 0.184 0.226 0.285 0.38 0.563
210.154 0.104 4.690.0164 1.64e-050.0004380.00181 0.00415 0.00752 0.012   0.0176  0.0244  0.0326  0.0424  0.0538  0.0673 0.0831 0.102  0.124  0.151 0.184 0.226 0.285 0.38 0.563
210.273 0.104 4.690.0164 1.64e-050.0004380.00181 0.00415 0.00752 0.012   0.0176  0.0244  0.0326  0.0424  0.0538  0.0673 0.0831 0.102  0.124  0.151 0.184 0.226 0.285 0.38 0.563
210.278 0.104 4.690.0164 1.64e-050.0004380.00181 0.00415 0.00752 0.012   0.0176  0.0244  0.0326  0.0424  0.0538  0.0673 0.0831 0.102  0.124  0.151 0.184 0.226 0.285 0.38 0.563
210.00010.104 4.690.0164 1.64e-050.0004380.00181 0.00415 0.00752 0.012   0.0176  0.0244  0.0326  0.0424  0.0538  0.0673 0.0831 0.102  0.124  0.151 0.184 0.226 0.285 0.38 0.563
220.296 0.112 4.690.0175 3.28e-050.0006940.00259 0.00562 0.00978 0.0151  0.0216  0.0294  0.0386  0.0494  0.0618  0.0763 0.0931 0.113  0.136  0.164 0.198 0.241 0.301 0.3960.577
220.04560.112 4.690.0175 3.28e-050.0006940.00259 0.00562 0.00978 0.0151  0.0216  0.0294  0.0386  0.0494  0.0618  0.0763 0.0931 0.113  0.136  0.164 0.198 0.241 0.301 0.3960.577
230.333 0.121 4.690.0187 6.23e-050.00106 0.00362 0.00746 0.0125  0.0188  0.0263  0.0351  0.0453  0.057   0.0705  0.086  0.104  0.124  0.149  0.177 0.212 0.256 0.317 0.4120.591
230.26  0.121 4.690.0187 6.23e-050.00106 0.00362 0.00746 0.0125  0.0188  0.0263  0.0351  0.0453  0.057   0.0705  0.086  0.104  0.124  0.149  0.177 0.212 0.256 0.317 0.4120.591
230.105 0.121 4.690.0187 6.23e-050.00106 0.00362 0.00746 0.0125  0.0188  0.0263  0.0351  0.0453  0.057   0.0705  0.086  0.104  0.124  0.149  0.177 0.212 0.256 0.317 0.4120.591
230.05270.121 4.690.0187 6.23e-050.00106 0.00362 0.00746 0.0125  0.0188  0.0263  0.0351  0.0453  0.057   0.0705  0.086  0.104  0.124  0.149  0.177 0.212 0.256 0.317 0.4120.591
230.231 0.121 4.690.0187 6.23e-050.00106 0.00362 0.00746 0.0125  0.0188  0.0263  0.0351  0.0453  0.057   0.0705  0.086  0.104  0.124  0.149  0.177 0.212 0.256 0.317 0.4120.591
230.05890.121 4.690.0187 6.23e-050.00106 0.00362 0.00746 0.0125  0.0188  0.0263  0.0351  0.0453  0.057   0.0705  0.086  0.104  0.124  0.149  0.177 0.212 0.256 0.317 0.4120.591
240.25  0.13  4.690.0199 0.0001130.00158 0.00496 0.00972 0.0158  0.023   0.0316  0.0414  0.0527  0.0655  0.08    0.0964 0.115  0.137  0.162  0.191 0.227 0.272 0.333 0.4280.605
240.04560.13  4.690.0199 0.0001130.00158 0.00496 0.00972 0.0158  0.023   0.0316  0.0414  0.0527  0.0655  0.08    0.0964 0.115  0.137  0.162  0.191 0.227 0.272 0.333 0.4280.605
240.179 0.13  4.690.0199 0.0001130.00158 0.00496 0.00972 0.0158  0.023   0.0316  0.0414  0.0527  0.0655  0.08    0.0964 0.115  0.137  0.162  0.191 0.227 0.272 0.333 0.4280.605
240.00010.13  4.690.0199 0.0001130.00158 0.00496 0.00972 0.0158  0.023   0.0316  0.0414  0.0527  0.0655  0.08    0.0964 0.115  0.137  0.162  0.191 0.227 0.272 0.333 0.4280.605
240.06260.13  4.690.0199 0.0001130.00158 0.00496 0.00972 0.0158  0.023   0.0316  0.0414  0.0527  0.0655  0.08    0.0964 0.115  0.137  0.162  0.191 0.227 0.272 0.333 0.4280.605
250.25  0.14  4.690.0211 0.0001950.00229 0.00664 0.0125  0.0196  0.0279  0.0376  0.0485  0.0608  0.0746  0.0902  0.108  0.127  0.15   0.176  0.206 0.243 0.289 0.35  0.4450.62 
260.05010.15  4.690.0224 0.0003250.00322 0.00872 0.0157  0.024   0.0335  0.0443  0.0563  0.0697  0.0846  0.101   0.12   0.14   0.164  0.191  0.222 0.259 0.306 0.367 0.4620.634
260.31  0.15  4.690.0224 0.0003250.00322 0.00872 0.0157  0.024   0.0335  0.0443  0.0563  0.0697  0.0846  0.101   0.12   0.14   0.164  0.191  0.222 0.259 0.306 0.367 0.4620.634
260.107 0.15  4.690.0224 0.0003250.00322 0.00872 0.0157  0.024   0.0335  0.0443  0.0563  0.0697  0.0846  0.101   0.12   0.14   0.164  0.191  0.222 0.259 0.306 0.367 0.4620.634
280.06910.172 4.690.025  0.00081 0.00599 0.0143  0.024   0.0349  0.0469  0.06    0.0743  0.0899  0.107   0.125   0.146  0.168  0.194  0.222  0.255 0.294 0.341 0.403 0.4970.664
280.22  0.172 4.690.025  0.00081 0.00599 0.0143  0.024   0.0349  0.0469  0.06    0.0743  0.0899  0.107   0.125   0.146  0.168  0.194  0.222  0.255 0.294 0.341 0.403 0.4970.664
280.22  0.172 4.690.025  0.00081 0.00599 0.0143  0.024   0.0349  0.0469  0.06    0.0743  0.0899  0.107   0.125   0.146  0.168  0.194  0.222  0.255 0.294 0.341 0.403 0.4970.664
280.235 0.172 4.690.025  0.00081 0.00599 0.0143  0.024   0.0349  0.0469  0.06    0.0743  0.0899  0.107   0.125   0.146  0.168  0.194  0.222  0.255 0.294 0.341 0.403 0.4970.664
300.1   0.197 4.690.0278 0.00178 0.0103  0.0222  0.035   0.0488  0.0634  0.079   0.0956  0.113   0.132   0.153   0.175  0.2    0.226  0.256  0.29  0.33  0.378 0.441 0.5330.694
300.154 0.197 4.690.0278 0.00178 0.0103  0.0222  0.035   0.0488  0.0634  0.079   0.0956  0.113   0.132   0.153   0.175  0.2    0.226  0.256  0.29  0.33  0.378 0.441 0.5330.694
320.111 0.224 4.690.0306 0.00353 0.0166  0.0327  0.049   0.0659  0.0833  0.101   0.12    0.14    0.161   0.184   0.208  0.234  0.262  0.293  0.328 0.369 0.418 0.48  0.5710.725
320.08580.224 4.690.0306 0.00353 0.0166  0.0327  0.049   0.0659  0.0833  0.101   0.12    0.14    0.161   0.184   0.208  0.234  0.262  0.293  0.328 0.369 0.418 0.48  0.5710.725
330.154 0.239 4.690.0319 0.00481 0.0206  0.039   0.0573  0.0757  0.0946  0.114   0.134   0.155   0.177   0.2     0.225  0.252  0.281  0.313  0.348 0.389 0.438 0.5   0.5890.74 
350.25  0.27  4.690.0346 0.0084  0.0306  0.0541  0.0763  0.0981  0.12    0.142   0.164   0.187   0.211   0.236   0.262  0.29   0.32   0.353  0.39  0.431 0.479 0.54  0.6280.77 
360.1   0.287 4.690.0359 0.0108  0.0366  0.0629  0.0872  0.111   0.134   0.157   0.18    0.204   0.229   0.255   0.282  0.311  0.341  0.374  0.411 0.452 0.501 0.561 0.6470.785
#modelpropFr$precision

#quantile_many=predict(modelpropFr, type = "quantile", at=c(.99))
#quantile_many

Al construir la figura para la regresión beta, una de las principales ventajas de utilizar este enfoque es que los cuartiles se calcula con una distribución beta. Por lo tanto, el margen de error NO baja de 0 y NO pasa de 1.

Evalua la siguiente figura en cada x hay una distribución beta, donde la linea roja representa ma mediana, las lineas verdes los cuartiles 25 y 75 y las lineas azules los cuartiles 5 y 95. NOTA que la distribución no es simétrica, y cambia a travez de la regresión.

library(ggplot2)
ggplot(dp2, aes(x=Number_of_Flowers, y=PropFR))+
  geom_point(width = 0.05, height = 0.01)+
  geom_line(aes(y=quantile_.05), linetype="twodash", colour="blue")+
  geom_line(aes(y=quantile_.25),linetype=2, colour="green")+
  geom_line(aes(y=quantile_.50), colour="red")+
  geom_line(aes(y=quantile_.75), linetype=2, colour="green")+
  geom_line(aes(y=quantile_.95), linetype="twodash", colour="blue")+
  ylab("Predicción de la proporción de frutos")+
  xlab("Números de Flores")+
  annotate("text", x=25, y=0.50, label="95th quartile", fontface="italic")+
  annotate("text", x=32, y=0.39, label="75th quartile", fontface="italic")+
  annotate("text", x=33, y=0.14, label="25th quartile", fontface="italic")+
  annotate("text", x=33, y=0.27, label="Median", fontface="italic")+
  annotate("text", x=35, y=-0.02, label="5th quartile", fontface="italic")+
theme(axis.title.y = element_text(colour="grey20",size=20,face="bold"),
        axis.text.x = element_text(colour="grey20",size=20,face="bold"),
        axis.text.y = element_text(colour="grey20",size=20,face="bold"), 
        axis.title.x = element_text(colour="grey20",size=20,face="bold"))+
  theme(legend.position="none")+
  rlt_theme
## Warning: Ignoring unknown parameters: width, height

ggsave("Graficos/Beta_number_flowers.png")
## Saving 7 x 5 in image

Aquí es una representación de las distribuciones beta en las x, número de flores. En rojo es para simular la distribución la proporción esperada de frutos en plantas que tienen 15 flores y la linea azul es para simular la distribución esperada de la proporción de frutos en plantas con 30 flores.
***

Evaluando la distribución de beta para valores específicos de x= Proporción de frutos por plantas basado en la cantidad de flores en la planta.

Seleccionar los diferentes valores de x y calcular el promedio y la varianza y convertir estos en \(\alpha\) y \(\beta\). Con estos parámetros se puede construir la densidad de la distribución por cada valor de x.

Se selecciona valores específicos para la visualizar la distribución, las plantas que tienen 25, 30 y 35 flores. Se tiene que reorganizar los datos para calcular el promedio y la varianza.

dpQ=dp2 %>% 
  select(c(1, 6:26))

dpQ2=dpQ%>% 
  filter(Number_of_Flowers== 13) %>% select(c(2:22))%>% t %>% as.data.frame
dpQ2$Quartiles=c(.01, 0.05, .1, .15, .2, .25, .30, .35, .4, .45, .5, .55, .6, .65,  .7,.75, .8, .85, .9, .95, .99 )

mean(dpQ2$V1, na.rm=FALSE)
## [1] 0.06453931
var(dpQ2$V1, na.rm=FALSE)
## [1] 0.01296751
#dpQ2


ggplot(dpQ2, aes(Quartiles, V1))+
 geom_line()

Usando la varianza calculada en el chunk anterior, se puede calcular el \(\alpha\) y el \(\beta\).

estBetaParams <- function(mu, var) {
  alpha <- ((1 - mu) / var - 1 / mu) * mu ^ 2
  beta <- alpha * (1 / mu - 1)
  return(params = list(alpha = alpha, beta = beta))
}
#mean(Smokers$Y2000P)
#var(Smokers$Y2000P)
estBetaParams(0.2757771,0.04226223)
## $alpha
## [1] 1.027498
## 
## $beta
## [1] 2.698331

Producción de los gráficos. Se observa que para las plantas que tienen 25 y 30 flores la densidad esta sesgada a la izquierda, en otra palabra la probabilidad de tener pocas frutos domina la distribuciones.

library(gridExtra)
a=ggplot(dipodium2, aes(PropFR))+
stat_function(aes(x = dipodium2$PropFR, y = ..y..), fun = dbeta, colour="red", n = 62,
                      args = list(shape1 =0.5418068, shape2 =3.135593))+
  annotate("text", label="x=25 flores", x=0.1, y=10)+
  ylab("Beta \nDensity")+
  xlab("Probabilidad de tener frutos")

b=ggplot(dipodium2, aes(PropFR))+
stat_function(aes(x = dipodium2$PropFR, y = ..y..), fun = dbeta, colour="blue", n = 62,
                      args = list(shape1 = 0.7549048, shape2 =2.949404))+
  annotate("text", label="x=30 flores", x=0.1, y=5)+
  ylab("Beta \nDensity")+
  xlab("Probabilidad de tener frutos")

c=ggplot(dipodium2, aes(PropFR))+
stat_function(aes(x = dipodium2$PropFR, y = ..y..), fun = dbeta, colour="black", n = 62,
                      args = list(shape1 = 1.027498, shape2 =2.698331))+
  annotate("text", label="x=35 flores", x=0.1, y=2)+
  ylab("Beta \nDensity")+
  xlab("Probabilidad de tener frutos")
tresDensidad=grid.arrange(a,b,c, ncol=1)

tresDensidad
## TableGrob (3 x 1) "arrange": 3 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (2-2,1-1) arrange gtable[layout]
## 3 3 (3-3,1-1) arrange gtable[layout]
ggsave("Graficos/tresDensidad.png")
## Saving 7 x 5 in image

“Activities reported in this website was supported by the National Institute of General Medical Sciences of the National Institutes of Health under Award Number R25GM121270. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.”